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We describe a Monte Carlo scheme which, in a single simulation, yields a measurement of the 
chemical potential of a crystalline solid. Within the isobaric ensemble, this immediately provides an 
estimate of the system free energy, with statistical uncertainties that are determined precisely and 
transparently. An extension to multiple occupancy ("cluster") solids permits the direct determina- 
tion of the cluster chemical potential and hence the equilibrium conditions. We apply the method 
to a model exhibiting cluster crystalline phases, where we find evidence for an infinite cascade of 
critical points terminating coexistence between crystals of differing site occupancies. 



The phase behaviour of crystalline materials is impor- 
tant in fields as diverse as solid state physics, soft mat- 
ter, mineralogy and pharmacology. For instance, metals 
and their alloys exhibit rich phase behaviour [1] (includ- 
ing novel features such as isostructural transitions [5]); 
colloids can self assemble into a variety of complex struc- 
tures with applications in photonics [3] ; many drug com- 
pounds exhibit crystalline polymorphism which can in- 
fluence their clinical function [4]. The staple simulation 
approach for predicting crystalline phase behaviour is via 
free energy estimates obtained by numerical integration 
along some path that connects the macrostate of inter- 
est to a reference state of known free energy E] . Such 
"Thermodynamic Integration" (TI) is popular because 
it is both conceptually simple and can be implemented 
with only a modest extension of the simulation frame- 
work needed for standard Monte Carlo (MC) sampling. 
However there are a number of respects in which it is less 
than ideal. The method hinges on the identification of 
a good path and reference macrostate. A 'good' path is 
short; but the reference macrostate (the choice of which is 
limited) may lie far from the physical macrostate of inter- 
est, entailing a large number of independent simulations 
to make the necessary link. A potentially more serious 
constraint on the path is that the derivative being mea- 
sured should vary slowly, smoothly and reversibly along 
it; if it does not the numerical quadrature may be com- 
promised. A phase transition en route is thus a particular 
hazard. Evidently one has to decide how many simula- 
tions are to be performed along the path and where. In 
so doing one must strike a suitable balance between min- 
imizing computation time and ensuring that no region of 
the path is neglected. This may necessitate a degree of 
trial and error. The uncertainties to be attached to TI 
estimates are also problematic. Use of simple numerical 
quadrature will result in errors. Error bounds have to 
aggregate the uncertainties (statistical and systematic) 
associated with different stages of the integration pro- 
cess. 



Beyond simple crystals, there is considerable interest 
in particles that self assemble via microphase separa- 
tion into periodically modulated nanostructures. Clas- 
sic examples include the lamellar and micellar crystalline 
phases encountered in surfactants and copolymers [7]. 
More recently, it has been discovered that when certain 
types of repulsive particles that lack a hard repulsive 
core are compressed to high density, multiple occupancy 
( "cluster" ) crystals are formed [8HTT] . Such coreless po- 
tentials serve as models for a wide range of soft matter 
systems such as star polymers, dendrimer and microgels 
in which particles can substantially overlap [T^]. To de- 
scribe cluster crystals one must allow for a crystalline 
lattice in which each lattice site can be occupied by mul- 
tiple particles. Let us suppose that such a crystal has iV c 
lattice sites, labeled i — l...N c and that site i is occu- 
pied by n c (i) particles (a "cluster"). Clusters are gener- 
ally bi- or poly-disperse, so the total particle number is 
N = n c{i) = N c n c , with n c the average occupancy. 

A particular problem for simulation is to determine the 
equilibrium values of n c , the lattice parameter a, the 
pressure P and the chemical potential p that correspond 
to some particle number density p — N/V and tempera- 
ture T of interest. As shown previously, measurement of 
the Helmholtz free energy F at fixed N c in the constant- 
(NVT) ensemble is insufficient in this regard [5]. Instead 
one has to estimate the lattice site or cluster chemical 
potential p c , given by N c p c = F + PV — pN, which van- 
ishes at equilibrium |13j . Doing so entails supplementing 
TI measurements of F with additional sampling of the 
chemical potential p (via the Widom insertion method) 
and the pressure P (via the virial) [5]. This process, or 
alternatively a direct estimation of the constrained free 
energy [TU] , then has to be repeated for a range of values 
of n c in order to pinpoint equilibrium. Accordingly it is 
cumbersome and laborious. 

Here we introduce a new MC simulation scheme that 
allows direct calculation of the chemical potential of crys- 
tals (and thence the free energy) from a single simulation. 



Extending the method to cluster crystals permits direct 
estimation of the cluster chemical potential, while his- 
togram reweighting techniques can be used to identify 
the equilibrium state without further simulation. Wc 
first describe the basic scheme for a simple crystal be- 
fore outlining its cluster solid-generalisation. 

The central idea is to construct (within the constant - 
NPT ensemble) a reversible sampling path between a 
lattice with N + M particles and another with N parti- 
cles. The relative probability of finding the simulation in 
the two states provides a measure of the Gibbs free en- 
ergy difference AG = pM. This yields the chemical po- 
tential p, from which the Helmholtz free energy density 
follows immediately as / = pp—P, with P the prescribed 
pressure and p the measured particle number density. 

To elaborate, consider the situation shown schemati- 
cally in Fig. [I] for a cubic lattice (though note the method 
is applicable to any Bravais lattice). A constant-TVPT 
ensemble Monte Carlo simulation [6j is to be found in one 
of two states a £ {0, 1}. For a = the system comprises 
a periodic box of volume V^ -* containing (m+1) xmx m 
unit cells of lattice parameter a. Each particle (a circle) 
is associated with a unique site of a fixed perfect lattice 
(a black dot): there are N particles in the cubic sub- 
volume of m 3 unit cells shown, and M — N/m particles 
in the remaining (rightmost) plane of unit cells. In the 
spirit of the phase switch method O US] , we write the 
position vector of each particle i in terms of the displace- 
ment Ui from its lattice site Rf \ i.e. r+ = + Ui, 
i = 1 . . . N + M. 

The switch to the a = 1 state comprises a reversible 
mapping in which the instantaneous particle displace- 
ments {ui} are re-associated with a second set of lat- 
tice sites such that r/ 1} = #f ) + u % . We set 
= Rf ] for i = 1...N. Thus for u = 1, the 
first N lattice sites form a periodic cubic system with 
volume = V^°'m/(m + 1), and the particles as- 
sociated with these sites retain their relative positions 
within the box under the switch. By contrast, for the 
remaining M particles, the change in environment under 
the switch is more radical; they leave the box altogether 
to become "ghost" particles, associated with fixed sites 
{^ (1) }, i = N + 1...N + M. Ghost particles are in- 
dependent (so the relative positions of the fixed sites is 
arbitrary) and they experience only a harmonic confin- 
ing potential (j> g (ui) whose amplitude is chosen to roughly 
match the average ghost particle displacement to that of 
real particles. 

Consider now the associated statistical mechanics. If 
we write the partition function of the system in state a 
as Z^ a \ then since ghost particles are independent, 
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FIG. 1: (Color On 



Schematic of the ghost particle 



switching scheme described in the text. 



ble) partition function of one ghost particle, with f3 — 
l/(k&T) as usual. The free energy change associated 
with the switch a = 1 — > follows as 



(3 AG = (3Mp = In 
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with p the chemical potential. In order to estimate AG, 
we supplement standard MC updates of the particle dis- 
placements and box volume with attempts to switch a. 
These are accepted with the standard (NPT) ensemble 
probability: p a = min[l, aexp (— /3(PAV + AE))], with 
a = y( 1 ~ cr )/y( <T ). In general, however, such switch at- 
tempts suffer low acceptance rates since the particle dis- 
placements {Hi} for the current a may not all be typical 
for the switched a. To deal with this we implement biased 
("umbrella") sampling [B], which enhances the probabil- 
ity of configurations (for each a) for which the instanta- 
neous switch cost k — f3(pAV + AE) is small. Specifi- 
cally we include a weight function r\ a (k) in the acceptance 
probabilities for all MC updates. Weights are obtainable 
via any of the standard techniques such as transition ma- 
trix or Wang Landau sampling. Their effects are unfolded 
from the sampling in the usual way [BJ at the end of the 
simulation. 

In this manner one accumulates the relative probabil- 
ity pW/p' ) = Z^ /Z^ ' of finding the simulation in the 
respective a states, from which the requisite chemical po- 
tential follows as 



where Z g = J duexjp(—f3(j3 g (u)) is the (exactly calcula- 



p = kBTlA-r 1 mfcW/pW) - \aZ g ] . (4) 
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Statistical errors in p are determined simply by the 
switching statistics and readily quantified via a block 
analysis. If the uncertainty in r = pW/pW is a, then 
that in p is 0(a/M). Since M is typically O(10 2 ), this 
bestows the method with high sensitivity. To test it wc 
have measured p for a Lcnnard- Jones fee crystal of sys- 
tem size m = 4. Interactions were truncated at r c = 2.9 
and a mean-field correction of the usual type [6J was ap- 
plied. As no previous estimates of chemical potentials for 
this system exist (to our knowledge), we use our results 
to calculate the absolute Helmholtz free energy density 
/ = fip — P, and compare with literature estimates deter- 
mined by TI. Only one such previous estimate quotes an 
associated uncertainty allowing for meaningful compar- 
ison [IB] . For the state point T = 2.0, p = 1.28, we 
find pP = 20.985(2). ftp. = 18.967(3), pf = 3.292(1), 
the latter comparing well with the estimate of Vega et 
al. Pf — 3.290(4), though our error bar is substantially 
smaller. 

We next address a challenging problem in the simula- 
tion of cluster solids. At some (T,p), equilibrium cor- 
responds to a particular value of the cluster occupancy 
n c and the lattice parameter a. But specifying p fixes 
neither of these parameters, e.g. in an fee cluster crystal, 
p = 4n c /a 3 , which can be realized by many combinations 
of n c and a. In a real system, N c changes to relax the sys- 
tem to equilibrium. However, in conventional simulation 
ensembles, the value of N c is constrained on accessible 
times scales by free energy barriers and does not fluc- 
tuate. As described above, one approach to determining 
equilibrium in these circumstances is to estimate the clus- 
ter chemical potential p c (n c ) via a laborious combination 
of TI, Widom particle insertion and virial sampling [5], 
while another is to directly estimate the constrained free 
energy via TI 10J. By contrast ghost particle switch- 
ing (or more precisely ghost cluster switching) provides a 
simpler, more efficient and elegant solution to this prob- 
lem. 

In seeking to apply the method, it is expedient to em- 
ploy an ensemble in which both n c and a are free to 
fluctuate-the constant p, P, T ensemble. This ensemble is 
rarely utilized in simulations because the extensive scal- 
ing of the entropy means that the partition function is 
finite for high pressures, diverges on approach to equilib- 
rium, e.g. Z(p,,P,T) ~ (P - Pcq) -1 as P -> P+, and is 
infinite for all lower pressures. This is no longer a prob- 
lem when we have a constraint of fixed N c . The partition 
function is then 

Z(p J ,P,T,N c ) = Y^ [ dVdEe siN > v ' E ^-^ E+pv -^ 

N 

The entropy will be extensive for large N c , 
S(N,V,E,N C ) = N c s(n c ,v c ,e c ) with n c = N/N c , 
v c = V/N c , e c = E/N c . The dominant contribution to 
Z comes from the maximum of the integrand, where 
—ftp = ds/dn c , PP = ds/dv c and P = ds/de c . Denoting 



these saddle point values with an asterisk, the extensive 
contribution to In Z is 

In Z(p, P, T, N c ) = N c [s(<,<,<) ~ P{< ~ iK + P<)] 

(5) 

In general, this is a non-equilibrium partition function 
because at equilibrium any two of (p, P, T) determine the 
third. To find the equilibrium condition, assume (p, P, T) 
is a set of equilibrium parameters, then so must n*, v* 
and e* be. This means that they are obtained by max- 
imizing the entropy S(N,V, E, N c ) over N c . From the 
extensive form of S above, and using that P = ds/de c 
etc, this shows directly that the combination in square 
brackets in ^ must vanish. 

The upshot of this analysis is that 
—kBTlnZ(p,P,T,N c )/N c vanishes at equilibrium, 
which is as expected given that by standard thermo- 
dynamic arguments this quantity can also be identified 
with the cluster chemical potential p c . On the other 
hand, Z(p, P,T, N c ) is also the weight of different N c 
values in a simulation where N c can fluctuate, and so the 
above equilibrium criterion tells us that at equilibrium 
these weights are to leading order independent of N c . 
Ghost cluster switching allows one to repeatedly add 
and remove a crystallographic plane of lattice sites, 
thereby circumventing the barriers between different 
values of N c . Measurements of the relative probabilities 
of macrostates with different iV c then directly probes p c . 

The implementation is technically similar to that de- 
scribed for simple crystals, except that the particle num- 
ber is permitted to fluctuate via insertions/deletions 
which for cluster solids are efficient owing to the lack of 
a repulsive hard core in the potential. For ghost sites, in 
addition to choosing the harmonic amplitude such that 
ghost particle displacements are similar to those of real 
particles, we impose a ghost chemical potential p g cho- 
sen to yield an average site occupancy close to that of 
real sites. The order parameter against which we bias 
to enhance the switch probability is extended from the 
single occupancy case to become k = P(PAV + AE ± 
(p — p g )Ng), with Ng the ins t ant aneous number of par- 
ticles associated with lattice sites i = N + 1 . . . N + Al. 
We sample the distribution of the observables (N, V, E) 
across the two values of a, from which we unfold the ef- 
fects of the biasing weights and the ghost particle free 
energy. This yields (inter alia), the volume distribution 
p(V\p,, P, T), which exhibits two peaks, one correspond- 
ing to cr = 1 and the other for a — 0. Equilibrium is 
signaled in a very simple fashion by the equality of the 
peak weights of p(V), as these have the same ratio as 
Z(p, P, T, N c = N) and Z{p, P,T,N C = N + M). 

To validate the methodology we have considered a 
prototype cluster solid: the generalized exponential 
model (GEM-4) whose interaction potential is u(r) — 
eexp((— r/c) 4 ). We determined P cq and p cq for T = 
1.1, p = 8.5 in the fee phase - a state point for which 



4 



prior TI data is available [5]. The initial simulation was 
performed for P = 114.5, /i = 29.7 and the resulting 
probability distribution p(V) was subsequently extrapo- 
lated in P and fi to yield both equal peak weights and 
a density matching the target value. Fig. [2] shows the 
resulting equilibrium form, obtained for P cq = 114.45(1), 
^ cq = 29.752(3). The associated cluster number is n c — 
17.470(5) and the fee lattice parameter a/ a = 2.018(1). 
We note that these results agree to within error with 
those of Mladek [T7] . Fig. [2] shows a portion of a snap- 
shot of an equilibrium configuration colored by cluster. 
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FIG. 2: (Color Online), (a) Equality of peak weights in p(V) 
that signifies equilibrium in the (p, P, T) ensemble (see text 
for the equilibrium parameters), (b) An equilibrium configu- 
ration viewed along the [100] direction. 

As an application we consider a recent proposal [11 J 
concerning the existence of a cascade of low tempera- 
ture critical points in the GEM-4 model. On the basis 
of ground state energy calculations and a phonon anal- 
ysis, these authors found a sequence of low temperature 
isostructural (fee) phase transitions on increasing density. 
At T = 0, each transition is such that the cluster number 
changes by unity from one integer value to the next, i.e. 
n c — > n c + 1. No theoretical evidence was found that 
any of these transition has an associated critical point 
at finite temperature, but the authors hypothesized that 
this should be the case. Indeed subsequent evidence for 
a critical point terminating the lowest density transition 
n c = 2 O 3 has been found using TI pH [18] . 

We have used the present method to search for fur- 



ther critical points in the model at higher density. We 
employed the known universal Ising form of the critical 
order parameter distribution to estimate the first four 
critical points [19j . The resulting critical point param- 
eters are listed in table [T] Surprisingly we find that T° 
is equal within error in each case as reflected in the in- 
dependence of the form of the density distributions at 
T = 0.04348 shown in Fig. [3} 



Transition 


p c P c p c T c 


2«3 
3«4 
4^5 
5^6 


1.239(4) 1.974(2) 2.9151(5) 0.0435(4) 
1.740(4) 3.879(2) 4.1878(4) 0.0435(4) 
2.257(5) 6.418(2) 5.4521(5) 0.0435(4) 
2.762(5) 9.575(3) 6.709(5) 0.0435(5) 



TABLE I: Estimated values of the critical density p c , pres- 
sure P c , chemical potential fi c and temperature T c , for the 
first four members of the infinite cascade of phase transitions 
having n c <-> n c + 1 at T = 0. 
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FIG. 3: (Color Online). Density distributions of the GEM- 
4 model corresponding to the near-critical point parameters 
listed in table [i] 

In summary we have introduced an efficient and ac- 
curate 'ghost particle switching' method for chemical 
potential determination in crystalline solids within the 
constant- NPT ensemble. The method, which circum- 
vents the need for integration to distant references states 
and its attendant pitfalls, requires only a single simula- 
tion at the state point of interest and yields statistical 
uncertainties directly and transparently. Such access to 
the chemical potential permits the direct determination 
of phase boundaries by matching of /i and P in the co- 
existing phases. An extension of the method to multiple 
occupancy crystals simplifies the problem of determin- 
ing their equilibrium parameters. As a demonstration 
of its power in this regard, we have studied the GEM-4 
cluster solid, uncovering the presence of a cascade of crit- 
ical points. More generally, the basic approach of ghost 
particle switching should be applicable to any system ex- 
hibiting periodic microphase separation such as lamellar 
or micellar crystals [7] , where the repeat unit can contain 
many individual particles. 
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